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I. INTRODUCTION 

Optical bistability is a manifestation of nonlinearity in 
optical systems which is well known in the laser physics 
community [T] [5] (see also [3 and [3] and references 
therein). It describes a situation in which there are two 
possible stable output light intensities for a single input 
intensity, and occurs when an optical medium with a non- 
Hnear refractive index is placed inside an optical cavity 
formed from two mirrors. The bistable behaviour results 
from the combination of the nonlinearity of the medium 
with the action of the feedback provided by the mirrors. 

A new addition to the family of systems displaying op- 
tical bistability has recently been demonstrated in ex- 
periments performed by the ultracold atom groups at 
Berkeley [S] and the ETH who found optical bista- 
bility in systems comprising of vapors of ultracold atoms 
trapped inside optical cavities which are driven by laser 
light. The atomic vapor acts as a dielectric medium and, 
despite being tenuous, can significantly perturb the light 
in a cavity with a small mode volume and high finesse if 
the cooperativity C^r is in the regime C^v = Ngf^l2K'y > 1, 
where N is the number of atoms, go is the single photon 
Rabi frequency, 27 is the atomic spontaneous emission 
rate in free space, and 2k is the cavity energy damping 
rate. The perturbation of the light by the atoms is non- 
linear and distorts the cavity lineshape away from being 
a lorentzian which is symmetric about the resonance fre- 
quency into one with an asymmetric shape. For large 
enough cooperativity the lineshape becomes folded over 
(see Fig. [2] below), so that for a certain range of frequen- 
cies there are three possible output light intensities (two 
stable, one unstable) from the cavity for a single input 
intensity. The experiments [5] and [5] exhibited this op- 
tical bistability as a hysteresis effect seen by chirping the 
laser frequency through the cavity resonance from above 
and below the resonance: a sudden jump in the inten- 
sity of light transmitted through the cavity was observed 



which occurred at two different frequencies, depending 
upon the direction of the chirp. 

An important difference between traditional laser sys- 
tems and the ultracold atom experiments [S] |S] is the 
origin of the nonlinearity. In the former case the nonlin- 
earity of the medium occurs in its polarization response, 
i.e. it arises from the internal degrees of freedom of the 
atoms. By contrast, in the ultracold atom experiments 
the detuning of the cavity from atomic resonance was 
large enough that the polarization response was in the 
linear regime. The nonlinearity was instead due to the re- 
sponse of the centre-of-mass wave function of the atoms: 
the atoms re-arrange their position distribution accord- 
ing to the balance between the dipole force applied by 
the intracavity light field (which forms a periodic lattice) 
and their zero-point energy. As a consequence, the depth 
of the optical lattice that forms inside the cavity in ex- 
periments like [S] and [B] is not fixed purely by the drive 
laser intensity, as is the case in standard optical lattices 
made by interfering laser beams in free space. Rather, 
when Cat > 1 the depth of the lattice is sensitive to the 
spatial distribution of the atoms trapped in the cavity, 
and, in turn, the atoms' center-of-mass wave function is 
sensitive to the lattice depth. This feedback nonlinearity, 
which leads to different amounts of transmitted/reflected 
light for a given input intensity depending on the spatial 
distribution of the atoms, has been previously employed 
to detect the presence of a single atom in a cavity [7 , as 
well as to monitor the motion of atoms trapped in cavi- 
ties [5]. More recently, there has also been considerable 
theoretical interest in the effect of the feedback nonlin- 
earity upon the many-body quantum state of ultracold 
atoms in cavities [5Hl3j. We note, in particular, the theo- 
retical work on self-organisation and related phenomena 
|14) . culminating in the experimental observation of the 
Dicke quantum phase transition |15j . 

Striking nonlinear phenomena also occur when ultra- 
cold atoms are trapped in standard "fixed" free-space lat- 
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tices [IS]. Of special interest to us here are the curious 
swallowtail loops that occur in the band structure (en- 
ergy versus quasi-momentum) of atomic Bose-Einstein 
condensates (BECs) in one-dimensional optical lattices. 
These loops have been studied theoretically by a num- 
ber of groups [I7H2TI in order to explain the breakdown 
in superfluidity observed in experiments where a BEG 
flows through a lattice [351 US] ■ The loops correspond to 
multiple solutions for the atomic wave function within a 
single band for a certain range of quasi-momenta. They 
can occur either around the boundaries of the Brillouin 
zone or the center, depending upon the band and the sign 
of the interactions. They manifest themselves physically 
via a dynamical instability that destroys the superflow. 
However, the nonlinearity responsible for the swallow- 
tail loops in the free-space lattices is provided by the 
interatomic interactions, which become important at the 
densities required for Bose-Einstein condensation. The 
loops occur when the strength of the interactions is above 
a critical value [T71 [TS], and therefore non-interacting 
atoms in an optical lattice do not display these instabili- 
ties. Our purpose in this paper is to investigate whether 
the cavity feedback nonlinearity associated with optical 
bistability in ultracold atoms can also lead to loops in the 
atomic band structure. As we shall see, the answer to this 
question is in general affirmative, and so band structure 
loops appear to be a robust phenomenon which appear 
whatever the source nonlinearity, although the structure 
and location of the loops does depend on the details of 
the nonlinearity. 

One consequence of loops in the atomic band structure 
is a hysteresis effect [H] if the quasi-momentum is swept 
through the band, and a consequent loss of adiabatic- 
ity even if the quasi-momentum is swept infinitely slowly 
[Mj . Recent experiments on a two-site lattice have con- 
firmed this scenario [25]. These effects have implications 
for experiments that use Bloch oscillations of cold atoms 
in optical lattices for high precision measurements, for ex- 
ample to determine the fine structure constant a |26j , to 
measure gravity [27l [28] , or to investigate Casimir-Polder 
forces [29] . The band structure hysteresis is reminiscent 
of the hysteresis described above in the context of the 
optical bistability experiments [S] and [S]. Indeed, as we 
shall show in this paper, for atoms in an optical cavity 
the two effects are different sides of the same coin, one 
being seen in the light and the other in the atoms. Of 
particular relevance, therefore, are two recent proposals 
[30l I3T] concerning Bloch oscillations of atoms inside op- 
tical cavities that rely upon the modification of the trans- 
mitted/reflected light arising from the feedback nonlin- 
earity as a non-destructive measurement of the Bloch os- 
cillations. The presence of loops will severely disrupt the 
Bloch oscillations. In the case where the nonlinearity 
has its origin in atom-atom interactions the loops can 
be eliminated by, e.g., using spin polarized fermions |29j 
or by reducing the interactions via a Feshbach resonance 
[32j . However, when the nonlinearity is due to the feed- 
back nonlinearity these methods do not apply, and one of 



our aims here is see if there are regimes where the feed- 
back nonlinearity is present and leads to a modification 
of the light, but remains below a critical value needed 
to form loops in the atomic band structure. In order to 
allow the origin of the nonlinearity to be clearly identi- 
fied, we shall only consider non-interacting atoms in this 
paper, but our calculations can be easily extended to in- 
clude interactions. 

The plan for the rest of the paper is as follows. In Sec- 
tion jll] we give a brief description of the system and in- 
troduce the mean-field hamiltonian describing cold atoms 
dispersively interacting with the single mode of a cavity. 
In Section jnij we derive a reduced hamiltonian describ- 
ing the nonlinear evolution for the atomic field after the 
adiabatic elimination of the light field. In Section |IV| we 
calculate the band structure by two different methods 
that solve for the spatially extended eigenstates (Bloch 
states) of the reduced hamiltonian, and show that the 
two methods give the same results. We find that for 
certain parameter regimes, the energy as a function of 
quasi-momentum develops loop structures. In Section [V] 
we recall the optical bistability in atom-cavity systems 
discussed by [5] [6] and make the connection between the 
loop dispersions and bistability. In Section|Vl]we develop 
an analytical estimate for the critical pump strength nec- 



essary to generate loops and in Section VII we test this 
result by illustrating how the band structure depends on 
laser detuning and laser intensity, i.e. the birth and death 
of loops as these parameters are varied. In Section [Vlllj 
we show that for quasi-momentum q the system can 
exhibit tristability (five solutions for the steady state cav- 
ity photon number). 

A convenient mathematical framework for describing 
bifurcations of a system whereby there is a sudden change 
in the number of solutions as a parameter is smoothly 
changed is catastrophe theory [33H37] . Catastrophe the- 
ory was applied to the problem of optical bistability in 
the late 1970s by a number of authors including Gilmore 
and Narducci [35] and Agarwal and Garmichael [? ]. In 
Section IX we relate the problem of atom-cavity multi- 



stability to catastrophe theory and show that the sys- 
tem under consideration can be described by swallow- 
tail catastrophes organized by an (unobserved) butter- 
fiy catastrophe and use this to understand multistability. 
This is followed in Section|X]by a discussion of the stabil- 
ity of these solutions and finally we conclude in Section 
XI There is also an appendix which summarizes some 
concepts of catastrophe theory. 



II. THE HAMILTONIAN 

Gonsider a gas of N ultracold atoms inside a Fabry- 
Perot optical cavity. The atoms interact quasi-resonantly 
with a single cavity mode of the electromagnetic field 
of frequency lUc, and it varies along the cavity axis as 
cos{ki,z), where = uJc/c This cavity mode is exter- 
nally pumped by a laser of frequency ujp. The relevant 
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frequency relations can be characterized with the two de- 
tunings 



Ac 



Wo 



Wo 



- Wc 
■ Wa 



(1) 

(2) 



where Ua is the atomic transition frequency. In the dis- 
persive regime, the occupation of the excited atomic state 
is vanishingly small and it can be adiabatically elim- 
inated. A one-dimensional hamiltonian for the atom- 
cavity system in the dispersive regime can then be writ- 
ten as [S HO] 

-SAcO^a -|- ihrj (a) — a) 



H 



dz ¥ 



2M dz 



a) acos^{kcz) 



(3) 



where and a{t) are the field operators for the 

atoms and the cavity photons, respectively. We have 
written the hamiltonian in a frame rotating with the 
pump laser frequency Wp, which leads to the appearance 
of the two detunings. The first term is just the free field 
evolution of the cavity mode. The second term represents 
the laser coherently pumping the cavity at rate 77, and the 
third term describes the atomic part of the hamiltonian. 
It consists of a kinetic energy part and a potential energy 
part. The potential energy term can either be understood 
as the atom moving in a periodic potential with ampli- 
tude {hgQ/ Aa){a^d), or, if combined with the first term 
in the hamiltonian, as a shift in the resonance frequency 
of the cavity due to the coupling between the atom and 
the field (see also Eq. ([5| below). 

The Heisenberg equations of motion for the atom and 
photon field operators can be obtained from the above 
hamiltonian. In this paper we are interested in proper- 
ties at a mean-field level, where operators are replaced by 
their expectation values and correlations between prod- 
ucts of operators are neglected. In other words, the 
light is assumed to be in a coherent state with amplitude 
a{t) = (a), and the atoms are assumed to all share the 
same single-particle wave function ii){x,t) = {"^)/Vn. 
From the Heisenberg equations these amplitudes obey 
the following coupled equations of motion [3D] 



da 
'dt 



+ ?7o»T-phCos2(a;) ) 



(4) 



iAc-iNUo / dxcos'^{x)\'ip{x)\'^ ~ k] a + r], 



(5) 



where we have scaled energies by the recoil energy — 
h'^kl/2M, and time by H/Er. The depth of the peri- 
odic potential seen by the atoms is then given by Uq riph, 
where riph = jap is the mean number of photons in the 
cavity, and is defined as 

C/o = hgl/[AaE^) . 



(6) 



The damping rate k (in units of E^) of the amplitude 
of the light field in the cavity has been added into the 



equation of motion in order to account for a markovian 
coupling between the cavity mode and a zero tempera- 
ture bath. We have also introduced the dimensionless 
length X = kcZ. In the above equations, any fluctuations 
induced by the reservoirs have been neglected. This is 
justified when considering relatively large quantum num- 
bers, for corrections see reference [T^ . 

In this paper we are interested in the band structure 
and this requires the quasi-momentum to be a good quan- 
tum number. Physically, this implies that we are study- 
ing delocalized atomic wave functions which cover many 
lattice sites (this is closer to the situation realized in the 
experiment [S] than that realized in [S] where the atoms 
were localized to just a few sites). We shall therefore 
expand the wave function ■0(x,t) in Bloch waves, as will 
be detailed in subsequent sections. It thus makes sense 
to normalize 'ip{x,t) over one period of the potential as 
/q \ip\'^dx = tt, and also evaluate the integrals appearing 
in the above equations over one period. In particular, the 
integral in Eq. ([s]) which determines the coupling between 
the atoms and the light will be defined by 

(cos2(a;)) = - r |^(a;)|2cos2(a;)da;. (7) 
Jo 



III. THE REDUCED HAMILTONIAN 

In this Section we shall eliminate the optical degrees 
of freedom from the hamiltonian in order to obtain a 
description only in terms of the atoms. This results in a 
nonlinear Schrodinger equation and an energy functional 
we shall refer to as the reduced hamiltonian. 

Setting da/dt = in Eq. ^ gives the steady state 
light amplitude in the cavity. 



a ■ 



V 



1 

A^~NUo{cos^{x)) 



(8) 



which leads to the following equation for the steady state 
photon number 

riph = ^ . (9) 

n^ + {Ac-NUa{cos-^ix))f 

In fact, this expression also holds to high accuracy in 
many time-dependent situations because k is typically far 
greater than any frequency associated with the evolution 
of the external degrees of freedom of ultracold atoms. 
The light field is then slaved to the atoms and "instan- 
taneously" adjusts itself to any change in ip{x,t). The 
steady state solution can then be substituted back into 
Eq. Q to give us a single, closed, nonlinear Schrodinger 
equation for the atomic wave function 
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The stationary solution ■0(a;, t) = ■ip{x)e of this 
equation leads to an expression for the eigenvalue of 
Eq. PI, 
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dx 



UoUph cos^(x)|'0(a;)|" 



(11) 



If the Schrodinger equation ( 10 ) were linear, then the 
eigenvalue /i would be the energy of the atoms in state 
tp. Furthermore, the functional (11) would serve as the 
energy functional from which this Schrodinger equation 
could be obtained as an equation of motion using Hamil- 
ton's equation 



ih 



dip _ SE[tP, ip* 
~dt ^ 



Sip* 



(12) 



i.e. E[ijj,ip*] — ^[ip,ip*] for a linear equation. However, 
this is not the case here. Instead, the eigenvalue /i is the 
chemical potential which is related to the true energy E 
via fi — OE/dN (a prominent example that illustrates 
this distinction is the Gross-Pitaevskii equation and its 
energy functional [H]). Using this fact, one can show 
that the true energy functional from which the equation 
of motion ( 10 1 can be derived is in fact JHl US] 



N 

^2 



dx 



arctan 
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as can be verified by applying Hamilton's equation (12 1. 
We shall refer to this functional as the reduced hamilto- 
nian. The first term represents the kinetic energy. The 
second term is an atom-light coupling dependent term 
that can be interpreted as follows. The phase shift of the 
steady state light field inside the cavity relative to the 
pump laser is, from Eq. ([8|, 
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This allows us to rewrite the reduced hamiltonian as 
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where /ph = tj'^/k is the incident photon current from 
the pump laser. Note that this hamiltonian looks similar 
in form to the effective quantum two-mode hamiltonians 
obtained in the optomechanical limit in [T^] and [13] , and 
also the ones describing bistability in [TU] . 



IV. BAND STRUCTURE 

From now on we specialize to Bloch wave solutions. 
We begin by describing two methods for calculating the 
Bloch waves and their energies. Agreement between the 
two methods will be demonstrated. 



A. Method 1: Energy extremization 

The first method, which adapts that detailed in [TS] for 
a static optical lattice, hinges on the observation that the 



potential in the Schrodinger Eq. ( 10 1 is periodic with pe- 
riod TT (in dimensionless units). Despite the nonlinearity, 
the Bloch theorem (31] applies so that the eigen- 

functions can be written as the product of a plane wave 
with wavenumber q, called the quasi-momentum, and an 
amplitude Uq{x) which is periodic with the period of the 
lattice 



Uq{x + Tr) = Ug{x) . 



(16) 



For the linear problem, the energies E{q) are arranged 
into bands separated by gaps. In the so-called reduced 
zone scheme for the band structure, q lies in the first 
Brillouin zone — 1 < g < 1, and the band energies are 
folded back into the first Brillouin zone. 

The periodicity of Uq{x) implies that it can be ex- 
panded as a Fourier series. The Bloch wave can therefore 
be written 



1Pg(x) 



A2nx 



(17) 



The notation for the n"^ expansion coefficient aq^n re- 
flects the fact that it depends on the chosen value of q. 
This expansion can now be substituted into the reduced 
hamiltonian Eq. (13), and the resulting function numer- 



ically extremized with respect to the parameters ag^„, 
maintaining the normalization of ipq{x) throughout the 
procedure. We take the parameters „ to be real for 
the same reasoning as given in |18| . The Fourier series is 
terminated at some n = R, which is determined by the 
convergence of the energy eigenvalues as R is varied. 

In Fig. [T] we plot the low lying part of energy spec- 
trum E[ipq] as a function of quasi-momentum resulting 
from the extremization. The values of k and N are very 
close to the values used in the experiment with rubid- 
ium atoms described in [B], and the other parameters 
are chosen so as to exhibit the interesting behavior to 
be discussed in the rest of the paper. The two panels of 
Fig.[l]differ only in the value of the pump-cavity detuning 
At;. Fig. [T]i shows a result familiar from linear problems 
involving quantum particles in periodic potentials, but 
Fig. [Tj^ shows a very different story: there are two other 
branches that together form a looped structure that is a 
manifestation of the nonlinearity of the reduced hamilto- 
nian. As will be discussed more below, the loop shown in 
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FIG. 1. Energy loops in the first band. The curves were obtained by extremizing the reduced hamiltonian ( |13[ ). For both (a) 
and (b) the laser pumping rj = 909.9 cjr, the cavity decay k = 350 wr, the atom-light coupling Uo = i^r, and number of atoms 
— 10*. In (a) the pump-cavity detuning was = 1350 a;R, which gives no loops, and in (b) it was Ac = 3140 wr which 
gives a loop symmetric about the band center as shown. As explained in the text, the number of photons in the cavity and 
hence the lattice depth change with q. For example, in (a) at g = we have riph = 0.06 and at g = 1 we have Wph = 0.68. In 
(b) at g = we have for the lowest branch Wph = 4.13, for the middle branch nph = 0.28, and for the upper branch nph ~ '2A. 
At g = 1 we have riph — 1.08. At the point where the middle and upper branches meet we have Uph ~ 0.58. 



Fig.[Ip belongs to the first band (in particular, it does not 
correspond to the second band because it only extends 
over part of the first Brillouin zone). Looped structures 
have been found before for BECs in static optical lat- 
tices [T71 [IHI . We will come back to the similarities and 



differences between our results and [T71 HH] in the next 
section. 

It is important to appreciate that, by virtue of the 
nonlinearity of the problem, the lattice depth riphC^ is 
different at each value of the quasi-momentum shown in 
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Fig. [T] Furthermore, the lattice depth is different for 
each of the curves even at the same values of q (except at 
degeneracies). In this sense, the band structure plots we 
display in this paper are non-standard because for static 
lattices the depth of the lattice is fixed for all values of 
q. In order to obtain the lattice depth at each point on 
a curve i?['0q] in Fig. fll one should take the wave func- 



tQe 



tion that extremizes the energy functional ( 13 1 at that 
point and enter it into the integral (cos^(a;))appearing 
on the righthand side of Eq. ^ for the photon number. 
This change in lattice depth with detuning is reported in 
Fig. [2] but for the reader's convenience we have included 
some values at particular points q in the caption of Fig.jl] 
The fact that, depending on the detuning Ac, there 
are either one or more steady state photon numbers for 
the cavity reminds us of the optical bistability observed 
in O |6] . There, as the cavity driving laser detuning was 
swept through the resonance, bistability was observed for 
certain pump strengths 77. The bistability derives from 
quantum effects [H] in the sense that it is due to changes 
in the atomic center-of-mass wave function. It is distinct 
from standard semi-classical optical bistability [3]. The 
connection between the loops in the atomic band struc- 
ture and optical bistability will be examined in detail in 
Section |Vj To complete this section we look at an alter- 
native method for determining the eigenfunctions of the 
effective hamiltonian which makes the connection with 
bistability more transparent. 



B. Method 2: Self-consistency equation 

In the second method, the strategy is to solve Eq. ([9| 
directly for the steady-state photon number. This equa- 
tion contains nph both explicitly on the left hand side 
and implicitly on the right hand side through the atomic 
wave function -05(0;, riph) appearing in the integrand of 
the integral 



(cos^(x)) 



co8^ {x)\il)q{x, nph)pda;. 



(18) 



The photon number is not a parameter set from the out- 
side (e.g. by the pump laser) but must be solved for 
self-consistently. In our notation for the wave function 
we have therefore included riph in the list of arguments 
rather than the list of parameters. The dependence of 
ipq(x,n-p\i) upon the number of photons is because the 
depth of the lattice in which the atoms sit is given by 
[/o'T-phj as can be seen directly from the Schrodinger 
equation Q. Therefore, Eq. ([9| must be solved self- 
consistently for nph, and we will often refer to it as the 
"self-consistency equation" . As mentioned in the intro- 
duction, the physical mechanism that gives rise to the 
feedback between the atoms and the light stems for the 
atom-light coupling, Eq. ^ [or Eq. ([S])]. The atoms' 
spatial distribution controls the phase shift suffered by 
the light on traversing the cavity, and hence the cavity 
resonance condition, and therefore the amplitude of the 



light in the cavity. At the same time, the amplitude of 
the light determines the depth of the lattice which influ- 
ences the atomic wave function. 

One class of solutions to the self-consistency problem 
is provided by the Mathieu functions [44j. In fact, these 
provide exact solutions because the Schrodinger equation 
Q for a particle in a sinusoidal potential is nothing but 
the Mathieu equation. Despite the fact that the ampli- 
tude of the sinusoidal potential in Q itself depends on 
the solution 'il)q{x,n-ph) of the equation, this amplitude 
evaluates to a constant because '0g(^;'^ph) appears un- 
der the integral given by Eq. ^ . All that is necessary is 
that the wave function that enters into the integral is the 
same as the wave function appearing in the rest of the 
Schrodinger equation. This is the case precisely when the 
self-consistency equation ([9| is satisfied. This leads us to 
a very important point: by virtue of the self-consistency 
equation (19]) , the photon number riph is completely equiv- 
alent, in the sense of the information it carries, to the 
wave function ijj. Said another way, the Mathieu func- 
tions are specified by only two quantities: the value of the 
quasi-momentum and the depth of the potential. Thus, 
once we set the quasi-momentum, which is a parameter, 
the wave function is entirely determined by f/o'^ph, where 
Uq is also a parameter. We shall frequently take advan- 
tage of the equivalence between tp and rip^ in the rest 
of this paper because it allows us to replace the wave 
function by a single number nph. 

We shall denote the Mathieu functions by Xm,q,npt,ix), 
where m is the band index. They satisfy Mathieu's equa- 
tion which in our problem takes the form 



^-hC/onphCOS (x) \Xni,q,npt,{x) 



=11. 



m,q,npiiXm,q,nph 



(19) 



Our notation for the Mathieu functions indicates the full 
parametric dependence with the exception of the atom- 
light coupling strength Uq. Note that in the Mathieu 
functions we list nph as a parameter, like q, because that 
is the role it plays in the Mathieu equation. We therefore 
have that 



(20) 



Mathieu's functions can, of course, be written in Bloch 
form. In this paper we focus on the first band and so we 
shall suppress the band index from now on. We therefore 
have Xq.riphi^) — 6xp(iga;)[/q_„pj^ (a;). Substituting this 
Bloch form into the time-dependent Schrodinger equa- 
tion Q as ^jjq{x,nph,t) = Xq,npt,ix)exp{-int) one ob- 
tains 

d \ ^ 

— Uq^npi,{x) + C/onphCos^(x)C/g,„pj^(x) 

= tJ-q,npi,Uq^npAx)- (21) 

This equation can either be solved numerically from 
scratch, or a package such as Mathematica can be used 
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which gives the Mathieu functions for each value of q, 
and J/oriph. For a particular choice of q, these Mathieu 
functions can now be used in ^ to find the value(s) of 
riph that give self-consistency. 

There are two main differences between method 1 and 
method 2 described above. Firstly, method 1 is a vari- 
ational calculation, whereas method 2 exploits the def- 
inition of the steady state photon number to obtain a 
single nonlinear integral equation ^ which must be sat- 
isfied by the atomic wave function. Secondly, in method 
2 we can explicitly set the band index whereas the vari- 
ational wave function used in method 1 is rather more 
general. In spite of these differences, we find that the 
two methods are in complete agreement (to within nu- 
merical accuracy) for all the parameter ranges we were 
able to test (for very deep lattices method 1 becomes un- 
workable because a very large number of terms in the 
Fourier expansion are required). We have also checked 
that the two methods agree for higher bands. 

It may at first seem surprising that two such seem- 
ingly different methods are equivalent. We emphasize 
that both stem from the non-linear Schrodinger equa- 
tion (10 1, which is just Mathieu's equation with a po- 
tential depth which is set self-consistently. Method 1 
minimizes an energy functional that is derived from this 
non-linear Schrodinger equation. In principle, one could 
minimize ansatze other than the Bloch functions we use 
here (e.g. localized functions), but these would not then 
satisfy the original non-linear Schrodinger equation ( 10 ) 
exactly. Method 2 is based upon the observation that 
wave functions that satisfy the self-consistency equation 
([9]) are precisely the required solutions of the non-linear 
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Schrodinger equation ( 10 1 providing we restrict ourselves 
to solutions which are Mathieu functions. Again, one 
could find other types of solution, but these would not 
satisfy Eq. pB. 



An interesting question to ask is whether the nonlin- 
earity of our problem mixes the linear bands, so that, for 
example the self-consistent first band is a superposition 
of Mathieu functions from different bands. This is not 
what we find. Instead, the solutions we have found all 
correspond to being from one or other band, but not a 
superposition. Method 2, in particular, allows us to ex- 
plicitly set the band index and we are therefore able to 
identify all three curves shown in Fig. [T] as belonging to 
the first band (we have also checked that the Mathieu 
functions corresponding to all three curves are orthog- 
onal to the Mathieu functions for the second band for 
the same three lattice depths). Actually, we do not find 
loops in higher bands for the parameter values consid- 
ered in Fig. [1] Although in this paper we restrict our 
attention to the first band, we have found that we can 
have loops in the higher bands as well. The calculation of 
energy dispersions using the self-consistent method is nu- 
merically less demanding and so we will continue to use 
the latter for the remainder of this paper. In the next 
section we discuss optical bistability and its relation to 
loops in the band structure. 



ri = 0.5 ri 
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FIG. 2. (Color online) The steady state photon number inside 
the cavity as a function of detuning Ac for the parameters k = 
350 Uq ~ i^a, A*' — 10*. Each curve is for a different value 
of the pump strength rj: thick blue rj = 0.5 77cr, dashed brown 
rj = 1.5 rjcT, thin magenta rj — 2.5 ria, and dash-dotted red 77 = 
3.5 77cr. As can be seen, as rj increases the lineshapes become 
more and more asymmetric and fold over at the critical pump 
strength ?7cr((? = 0) = 770 = 325 wr. The atomic wave function 
corresponds to the first band with g = 0. 



V. BISTABILITY AND LOOPS 

As mentioned in the Introduction, optical bistability 
was discovered in the ultracold atom experiments [5] and 
^Bj via a hysteresis effect as the detuning of the pump 
laser was swept from above and below the cavity reso- 
nance. This effect can be described theoretically by using 
the self-consistency equation to calculate the number of 
photons riph in the cavity at each value of the detuning 
(the number of cavity photons can be measured directly 
from the photon current transmitted through the cav- 
ity which is given by nphu). The results are plotted in 
Fig. [2] for different values of the pump strength 77. In 
the absence of atoms the cavity lineshape is a lorentzian 
centered at Ac = with a full width at half maximum 
2k. At small pump intensity, the presence of the atoms 
shifts the center of the resonance by NUo{cos'^{x)) while 
the shape is largely unaltered. But as r] is increased, the 
lineshape curve becomes asymmetric and eventually folds 
over when rj is above a critical value rjcriq = 0) = rjo, in- 
dicating multiple solutions and hence bistability. 7ycr('z) 
depends on the quasi-momentum and r]o is defined as its 
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FIG. 3. (Co lor online) The energy given by the reduced hamil- 
tonian ( 13 1 as a function of detuning Ac for the parameters 
q — 0, K = 350 ojR, Uo — io-R, N — W^. Each curve is for a 
different value of the pump strength rj: thick blue rj = 0.5 T^cr, 
dashed brown rj = 1.5 77cr, thin magenta rj = 2.5r;cr, and 
dash-dotted red r) = 3.5 rjcr- The critical pump strength is 
Tjo = 325 uJTi- For r/ > rjo swallowtail loops develop corre- 
sponding to bistability. The loops grow in size as rj increases. 



value at q = 0. 

In the folded over region, only one the solutions (corre- 
sponding either to the highest or the lowest photon num- 
ber) is seen at a time, depending upon the direction of 
the sweep. The middle branch cannot be accessed using 
this experimental protocol and is in any case unstable, 
as will be discussed at more length in Section pC} 

Note that in Fig. [2] we have chosen the quantum state 
of the atoms to be the q ^ Bloch state. In fact, this is 
the case for all the plots in this section because that is 
closest to the situation in the experiments that have been 
conducted so far. One of the main points of this paper is 
essentially to examine the extra degree of freedom con- 
ferred by the quasi-momentum. For atoms in ordinary 
optical lattices the quasi-momentum can be controlled 
by accelerating the lattice (rather than the atoms) via a 
phase chirp [iS] . An accelerated lattice is hard to achieve 
inside a Fabry-Perot cavity, but if the atoms have a mag- 
netic moment one can instead subject them to an inho- 
mogeneous magnetic field so that a force is applied (or, 
of course, in a vertically oriented cavity the atoms will be 
accelerated by gravity) . As is well known from the theory 
of Bloch oscillations, under a constant force F the quasi- 
momentum evolves according to the Bloch acceleration 




FIG. 4. (Color online) The double well structure of the 
energy of the reduced hamiltonian ( 13 1 as a function of 
cavity photon number nph. Each curve is for a differ- 
ent value of the detuning Ac. The values are Ac = 
1600, 2000, 2400, 2800, 3200, 3600 tJR. The arrow indicates 
how the curves evolve as Ac increases. The top most curve 
(blue), and the bottom most curve (yellow), have only one 
minimum, whereas the rest of the curves have two minima 
(the inset shows a zoom-in of the curves close to riph = 0) 
indicating bistability. Consequently, bistability only occurs 
for a certain limited range of Ac. The other parameters are 
q = 0, K, — 350 OJR, Uo — cjR, 
rjo = 325 a;R. 



A'^ = 10 , r; = 3.5 ?7o, where 



theorem 



lit) = qo 



Ft 



(22) 



where q and t are expressed in the dimensionless units 
given in Section [Ilj and go is the quasi-momentum at 
time t — Q. The Bloch acceleration theorem requires 
that the evolution be adiabatic, and the implications of 
this for intra-cavity optical lattices have been discussed 
in [3T], albeit without loops. The effect of a constant 
force is thus to sweep the system through the band at 
a constant rate and, in principle, any quasi-momentum 
can be achieved by switching off the magnetic field after 
an appropriate time delay. 

Figure [3] depicts the energy versus detuning curves cor- 
responding to the photon number versus detuning curves 
shown in Fig. [2] In the bistable regime, the energy curves 
in Fig. [3] develop swallowtail loops. This can be under- 
stood in terms of the familiar connection between bista- 
bility, hysteresis, and the change in the energy manifold 
described in detail by [T^. Consider one of the curves 
in Fig. [3] where there is bistability, e.g. the curve with 
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r] = 3r]Q. For Ac values to the left and right of the swal- 



lowtail loop, the energy functional Eq. ( 13 1 has a single 



extremum corresponding to a particular wave function 
ipq^XjUpii). In the bistable region, the energy functional 
has the structure of a double-well, furnishing three ex- 
trema: two minima and one maximum, that give the 
three branches of the loop corresponding to three differ- 
ent wave functions. This double-well structure is shown 
in Fig. I4] as a function of nph for different values of de- 
tuning Ac- Note that, as already observed above, the 
self-consistency equation ^ provides a direct mapping 
between the wave function and riph. 

Figures[2][3] and|4]all show different aspects of hystere- 
sis as Ac is swept either from above or below the cavity 
resonance. It is enlightening to see how it arises in Fig.|4] 
If the detuning is swept from below the resonance then 
initially there is a single solution for riph, given by the 
minimum in the reduced energy functional which occurs 
at the very left hand side of Fig. |4] as best seen in the 
inset. When Ac is increased another solution appears at 
a larger value of nph- However, this state of the system 
is not realized (for this direction of the detuning sweep), 
even when it becomes the global minimum, until the en- 
ergy barrier between the two solutions vanishes and the 
left hand minimum disappears. The system then jumps 
to the new minimum at a larger value of Tiph. The reverse 
happens when Ac is swept in the other direction. This 
hysteretic behavior is corroborated by Figs [2] and [Sj 

In the Section VI a method for determining rjcriq) is 



described. Generally, this requires a numerical computa- 
tion, but for small values of the intracavity lattice depth 
an analytical expression can be worked out. It turns 
out that the dependence of r]cT on the atomic state quasi- 
momentum can be used to explain the loops in the energy 
quasi-momentum plots (Fig. [T]). 



VI. CRITICAL PUMP STRENGTH FOR 
BISTABILITY 

A. Conditions for bistability 

Returning to the cavity lineshape shown in Fig. [2j we 
recall that as the pump strength 77 is increased the steady 
state photon number in the cavity can exhibit bistability 
for a certain range of detuning Ac. Bistability first devel- 
ops at a single value of the detuning, which we denote by 
Ac = Aq. The critical pump strength at which this bista- 
bility at Aq occurs is r]cr{q), and in this section we want 
to calculate it. Let us first re-write the self-consistency 
equation (23) as 



"-ph 



1 



1 



-NUof{U„7l^h,q) 



(23) 



where rimax = '7^/'*^ is the maximum number of the pho- 
tons that can be in the cavity at steady state. In order to 
reinforce the idea that the wave function and the number 
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FIG. 5. (Color online) Plot of the atom-light overlap integral 
/(nph,g) = (cos^(a;)), first defined in Eq. ([t]), as a function 
of the cavity photon number riph. The atomic wave function 
is taken to be the q = Bloch wave of the first band, and 
the atom-light interaction is set at Uo ~ 5cjr. Note that the 
maximum value /(nph, q) can take is one half, irrespective of 
the values of Uo and q. As riph — > 00 we find that /(wph, q) — >■ 
0. 



of photons are really equivalent quantities, we have re- 



placed the notation for the integral 
in pSl) by 



'cos (a;)) appearing 



fiUonph,q) 



(24) 



This function is plotted in Fig. [5] for blue detuning 
(Aa > 0) where we see that as the lattice gets deeper the 
atomic wave function adjusts to become more localized 
on the lattice nodes and reduce the overlap between the 
light and the atoms. Furthermore, the steep gradient at 
shallow lattices implies that the system is more sensitive, 
and so more nonlinear, at small photon numbers. 



It is instructive to solve Eq. ( 23 1 graphically as the 
intersection between two functions of Tiph, as shown in 
Fig. [6] The left hand side is a straight line whose gradient 
is 1/rimax- For very small nmax the gradient is very large 
and there is only one solution close to the origin. As rimax 
is increased the gradient is reduced and the straight line 
just grazes the curve at a critical value of n,„ax at which 
there are now two solutions. Increasing nmax further, 
there is then a range of values of nmax for which there 
are three solutions. Finally, for large nmax there is only 
one solution again. When three solutions exist for certain 
values of rimax, the system becomes bistable and a plot 
of the input intensity (proportional to rimax) versus the 
output intensity (proportional to riph) has the classic s- 
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FIG. 6. (Color online) A graphical solution of the self- 
consistency equation in the form (23 1. The blue curve rep- 
resents the right hand side of Eq. ( 23 1 for typical values of 
the cavity parameters (g — 0). The red dash-dotted, green 
dashed and black straight lines represent the left hand side of 
the equation plotted for different values of rimax; they inter- 
sect the blue curve at one, three, and one points, respectively. 
The blue curve tends to a finite value at Tiph — which is set 
by the fact that for nph = we have f = 1/2. 



shaped form shown in Fig. [7] This picture suggests a 
convenient way to determine the conditions for bistabihty 
because the two points where the curve turns over delimit 
the bistable region. These points satisfy 9?T-max/5nph = 
0, giving 



+ (Ae - NUoff 

- 2nph (Ae - NUof) NUo 



Of 



0. 



(25) 



This equation can be solved numerically for Uph for dif- 
ferent values of A^, assuming that k, N, and Uq are fixed. 
As expected from Fig. [t] (see also Fig. 14 ) , depending on 
the value of Ac, there are either zero, one, or two values 
of riph that satisfy Eq. ( [25] ). For large values of Ac there 
are no solutions. As Ac is reduced, a single solution for 
Uph suddenly appears at Ac = Aq, which, by substitut- 
ing this value of Uph into the self-consistency equation 
(231, gives us rjcriq). Reducing Ac further, this single so- 



lution immediately branches into two solutions for nph. 
Referring to Fig. [7j we see that these two solutions for 
Uph define a range of values of n-max, and hence also of 
r], where bistabihty occurs. We can find this range of 
values of rj by inserting the two solutions for np^ into 
the self-consistency equation to give us rji and 772 ■ When 




max 



FIG. 7. (Color online) Input intensity vs output intensity for 
a bistable cavity system. In this example the atomic wave 
function in the cavity is in the q = state and k — SSOujb., 
Ac — ISOOojR, and Uo ~ (^r. The points where the curve 



folds over are given by the solution of Eq. ( 25 1 



r?i < 77 < ?72 bistabihty occurs. Note that the values of 
Vcriq), and 772, all depend on the state of the atomic 
wave function and hence are different for different quasi- 
momenta. This dependence on quasi-momentum is what 
lies behind the existence of loops in the band structure. 



B. Critical pump strength in shallow lattices 



In the regime of shallow lattices, the bistabihty condi- 



tion given by Eq. ( 25 1 can be solved analytically. First 
consider the critical pump strength of the q = case, for 
which we use the notation rici {q = 0) = 770 • As described 
in 6 , and as can be seen from Fig. [5] for small lattice 
depths we can linearize the atom-light overlap integral as 



fiUonph,q 0) 



Upnph 
16 ■ 



(26) 



Substituting this into the self-consistency equation (|9|, 
we obtain a cubic equation in ?T,ph, which is reminiscent 
of the classical Kerr nonlinearity in a medium with an 
intensity dependent refractive index |H (note that when 
we come to numerically solve Eq. ( |25[ ) below in this sec- 
tion and in the rest of the paper, we are going beyond the 
classical Kerr effect). The condition (25 1 for bistabihty 



in this limit then reduces to the solution of a quadratic 



11 



equation in riph, 

Snli^N^U^ - 32nph Arc/2 (iVC/o - 2 A,) 

+ 64((7VL/o-2Ac)2+4k2) =0. (27) 

The vanishing of the discriminant of the above equation 

and 



requires that Aq = NUq/2 - V^k, tjq = ^ ^^^ui ' 

no = ^^4^^ • In the last expression, no is the numb 
photons in the cavity at the critical point Ac = Aq. 



C. Critical pump strength as a function of 
quasi- momentum 

We now extend the above analysis for shallow lattices 
to include non-zero quasi-momentum. Expanding the 
atomic Bloch state in a Fourier series, and assuming shal- 
low lattice depths, we can truncate the series after three 
terms 

^,(x,7iph,t) = e''^^(co(t) + ci(0e*'^+C2(i)e-''2-). (28) 

In this state one can explicitly calculate {cos'^{x)) as 
1 1 



(cos (x)) = - + -(RecoC2 -\-Recocl) 



--(X{t)+Y{t)). 



(29) 



Rewriting the equations of motion Q and ([5| for the 
newly defined variables one finds 



d^F 
dt 



+ {Aq + AfX + {q + l)[/oa*a|coP = 0, 



2 + (4g - AfY ~{q- l)C/oa*«|coP = 0, 



-=^^A^-^—-^ — iX + Y)-nja + V■ 

The atomic state has therefore been mapped onto two 
coupled oscillators X and Y. The oscillators are not cou- 
pled to each other directly, but do interact through the 
light field a which acts as a driving term for both of 
them. The above equations resemble Eqns (3) and (4) 
of [6], and introduce an analogy to optomechanics [46] , 
Solving these equations at steady state gives 



{Aq + AfX+- 



(4g-4)2y-- 



K-i{A^- NUo/2{X + Y)) 
iq + l)UoV^ 

Ac - NUo/2iX + Y) 

{q-l)UoV^ 



\co\ 



0, 



k2 -f Ac - NUo/2{X + Y) 



where Ac = Ac — NUq/2. Combining the steady state 
solutions into a single equation for the variable p = X + 
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FIG. 8. (Color online) Comparison between exact numerical 
calculation (dots) and analytical estimate (line) for the critical 
pump strength rjcr at which loops appear as a function of the 
quasi-momentum. The values of the parameters are Uo = ujr, 
N = 10*, and k — 350a;_R. The analytical estimate is from 



Eq. (31 1 which is accurate for small lattice depths. Note that 
the agreement is good for quasi-momentum close to g = 0. 



Y , and assuming jcop 
P 



1 gives 



1 



I 1^ A, _ NUa _ NUop V 
\ K, 2k 2k J 



(30) 



Comparing Eqns (30) and (23) 



TI ^ 

where fimax = siq^^-r 
in the limit when f — 1/2 ~ C/onph/16, we hnally obtain 
an expression for the critical pump strength above which 
bistability occurs as a function of q 



/8K3(l-g2) 



(31) 



where we remind the reader that the frequencies in this 
expression are in units of the recoil frequency un. This 
estimate for rjciio) is compared to the full numerical so- 
lution of Eq. (251 in Fig. [sj The parameters are such 
that the maximum intracavity depth is only of the order 
of one atomic recoil energy E^, , and hence the approxi- 
mation agrees well with the numerical calculation. The 
agreement deteriorates as q 1. This is due to the fact 
that the above two mode approximation fails at q = 1 



because the coefhcient cq in Eq. (28) is equal to zero at 
the edge of the Brillouin zone. 

Let us now connect the above results to the phe- 
nomenon of loops in the band structure. Because Tycr, 
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and also ryi, and 772 depend on the value of g, as we vary 
q we expect that the conditions required to have bista- 
bility won't necessarily be met over the entire Brillouin 
zone. That is, as q is varied, rj may no longer lie in the 
range 771(g) < rj < 772(g)- In that case, we expect any 
additional solutions to form closed loops extending only 
over part of the Brillouin zone, rather than entire bands 
covering the whole Brillouin zone. The dependence of 
77cr upon q given by Eq. (31) suggests that for shallow 



lattices the loops will form first at the edge of the Bril- 
louin zone and then propagate inwards as 77 is increased. 
Looking back at Fig. [T] which shows the energy plotted 
as a function of quasi-momentum at a fixed value of 77 
and Ac, we see a loop centered at g = 0, but which does 
not extend out to g = 1, in apparent contradiction to 



what is predicted by Eq. (31 ). This is because the lattice 
in Fig. [T]is too deep for Eq. (31 1 to apply. In the next 
section we shall examine how loops appear and disappear 
and, in particular, we will see that loops are indeed born 
at the edges of the Brillouin zone. 



VII. THE BIRTH AND DEATH OF LOOPS 

In this Section we examine how loops appear and dis- 
appear in the band structure as the detuning and pump- 
ing are varied. We have already seen in Sections |IV| and 
[V] how multiple solutions and swallowtail loops develop 
as the detuning from the cavity resonance is varied, but 
this was for fixed quasi-momentum g. Here we include 
the quasi-momentum dependence. In Fig. [9] we plot the 
evolution of the loop structures that appear in the band 
structure as Ac is varied. The detuning increases from 
the top to the bottom panel. In the plots, the pump 
strength is fixed at 7; = 2.8770 and we see that for small 
detunings, when bistability initially sets in, swallowtail 
shaped loops appear at the outer edges of the Brillouin 
zone. As the detuning is increased the swallowtail loops 
from the two edges move closer and merge. Initially, the 
merged loop lies partly below the lowest band, but as the 
detuning is further increased it moves up and separates 
from the lowest band. The loop then shrinks in size and 
vanishes. 

One important point to notice is that the swallowtail 
is qualitatively different from 
T8] for an interacting BEG in 



loop in plot (b) of Fi g. [9 
the ones obtained in [TtT 
an optical lattice because in our case the energy disper- 
sion continues to have zero slope at the band edge even 
when the loops have been formed. The nonlinearity in 
|17|, I18j . which is due to interatomic interactions, has 
quite a different form to that considered here. For exam- 
ple, the nonlinearity arising from interactions is spatially 
dependent due to the variation in density of a BEG in an 
optical lattice, whereas the nonlinearity considered here 
does not have a spatial dependence because it appears 
under an integral over space, see Eq. ([7|. 

In Fig.[TO]we plot the same thing as Fig. [9] except that 
we have reversed the sign of the atom-light coupling Uq 
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FIG. 9. (Color online) The birth and death of band struc- 
ture loops as the laser-cavity detuning Ac is varied, for the 
case when the laser is blue-detuned from atomic resonance 
(An > 0). Ac increases as one goes from (a) to (e) as follows: 
1500 UR, 2100 UR, 2600 UR, 3100 ur and 3600 a;R. The rest of 
the parameters are k = 350 ojr, Uo = ujr, N = 10*, 77 = 2.8 770 
and rio — 325 ojr. 
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FIG. 10. (Color online) The birth and death of band struc- 
ture loops as the laser-cavity detuning Ac is varied, for the 
case when the laser is red-detuned from atomic resonance 
(Aa < 0). Ac increases as one goes from (a) to (e) as follows: 
-8500 WR, -7900 tjR, -7400 WR, -6900 a;R and -6400 wr. 
The rest of the parameters are k = 350 wr, Uo = —cun, 
N = 10*, r; = 2.8770 and rjo = 325 a;R. 



SO that it is negative. Experimentally, this is the case 
when Aa < 0, i.e. the pump laser is red detuned from 
the atomic resonance. Note that the effect of the sign 
flip upon the potential term J/o'T.ph cos^(a:) occuring in 
the atomic Schrodinger equation is equivalent to a spa- 
tial translation of tt/2. This transforms the atom-light 
overlap integral ^ as 



(cos^(a;)) ^ 1 - (cos2(.t)) 



(32) 



where the left hand side refers to the Uq < case, and the 
right hand side to the Uq > case. The self-consistency 
equation therefore becomes 



77- 



k2 + (A, + N\Uo\ - N\Uo\f{\Uo\npi„ q)Y 



(33) 



Figure [TT] plots the evolution of the band structure as 
the external laser pumping 77 is varied, for fixed detuning 
Ac. We see that as 77 increased the loops first appear as 
swallowtails at the edges of the Brillouin zone (see inset 
in Fig. |11^ ), in agreement with the predictions of Eq. 

dsTl). 



Figures [12] and [I3| are both 3D plots of the loops as 
functions of A^ and 5, but each one covers a different 
range of Ac. In the first plot the merging of the swal- 
lowtail loops into the band center loops is shown. In the 
second plot, as Ac increases the band center loops move 
upwards in energy and separate from the ground band, 
shrink in size and eventually disappear. 

As will be demonstrated in the next section, in cer- 
tain parameter regimes the spectrum may qualitatively 
change compared to the above. In particular, we show 
that for some g ^ 0, and for sufficiently larger pump 
strength 77, wc can achieve tristability. Moreover, we 
show how this multistable behaviour can be understood 
in terms of catastrophe theory. 



VIII. TRISTABILITY 

Thus far we have seen that there are regions of the 
parameter space {Ac, 77, q, Uq} where the self-consistency 
equation ^ applied to the first band admits either one 
or three solutions. In the latter case we have bistability. 
It is natural to ask whether there are other regions of 
parameter space where even more simultaneous solutions 
can occur? A recent paper discussing two-component 
BECs in a cavity has found regions of parameter space 
that support tristability TT and in this section we want 
to examine whether tristability is possible in the ordinary 
one-component case but with finite values of the quasi- 
momentum. 



In Fig. 14 we show the region of {Ac, 77} space where 
bistability occurs. This plot is for fixed values of Uq and 
q. In particular, for q ^ Q we at most find bistabil- 
ity. The crescent shape in Fig. [M] demarks the region 
with three solutions: the number of solutions changes 
by two as one crosses its boundary. The location of 779, 
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FIG. 11. (Color online) The birth and death of band structure loops as the pumping rate r] is varied. In this figure the detuning 
is held constant at Ac — 2900 ujr. The value of rj increases from (a) to (e) as follows: O.Sryo, 2r7o, 3?7o, 4?7o, 5r;o, where rjo = 325 ojh 
as usual. The inset shows a zoom-in for rj — 0.5 770, illustrating that as rj is increased, the loops are born at the edges of the 
Brillouin zone. 
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FIG. 12. (Color online) Energy as a function of quasi- 
momentum q and detuning Ac. At smaller values of Ac, the 
swallowtail loops occur in pairs close to the edges of the Bril- 
louin zone at q = ±1, and as the detuning is increased they 
propagate inwards and merge as shown in this plot. Ac in- 
creases out of the page. Parameters are k = 350 cjr, Uq = 
N = 10^, 7? = 2.8 ?7o, and r?o = 325 cjr. 
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FIG. 13. (Color online) Energy as a function of quasi- 
momentum q and detuning Ac. Ac increases out of the page. 
For larger values of Ac, the band center loops move up in en- 
ergy and do not touch the lower band (shaded red in the plot). 
Eventually, for still higher values of detuning, they shrink and 
disappear as shown in this plot. Parameters are k = 350 ojr, 
Uo = a;R, N = 10*, r) = 2.8770, and t)o = 325 a;R. 



the smallest pump strength for which bistability occurs 
when g = 0, is indicated by an arrow on the vertical axis 
in order to make contact with the discussion of Sec. fVl 
However, when we allow non-zero values of q we find that 
we can indeed have regions of tristability, due to the pres- 
ence of five solutions, which occur inside the swallowtail 
shaped region in Fig. [151 which is plotted for q = 0.95. 
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FIG. 14. (Color online) Bifurcation structure of the solutions 
to the self-consistency equation Eq. ([ojl in the {77, Ac} plane 
with q — Q, Uo — a;R, k = 350clIr,A'' = 10*. The numbers 
on the plot indicate the number of solutions that exist for the 
steady state photon number in the cavity. The critical value 
of the pumping t^o = '7cr(? = 0) for bistability for these pa- 
rameters is indicated by the arrow. Inside the crescent shaped 
region the system supports three solutions (one unstable), i.e. 
it is bistable. 



Fig. [16] plots the corresponding photon number versus 
laser pumping curve for a fixed value of Ac, and illus- 
trates how, as the laser pumping strength is changed, 
the system goes from one solution, to three solutions, to 
five solutions, back to three solutions and finally back to 
one solution again. This curve is calculated for a verti- 



cal slice through the parameter space shown in Fig. 15 
We give an example of the band structure when there is 
tristability in Fig. [iTj 

So far we have conducted a rather ad hoc explo- 
ration of the four-dimensional parameter space given by 
{Ac, 77, (7, Uq\. Furthermore, in the two-dimensional slices 
shown in Figures [T4| and [TS] we have glimpsed snap shots 
of a rather complex looking structure of solutions. We 
are therefore led to ask whether there is any order in this 
complexity? Are the geometric structures seen in the 
plots random, or is there in fact an underlying structure 
that organizes them? In order to make progress with un- 
derstanding multistability it would be useful to have a 
more systematic framework for analyzing our solutions 
and just such a framework is provided by catastrophe 
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FIG. 15. (Color online) Bifurcation structure of the solutions 
to the self-consistency equation Eq. ^ in the {r/, Ac} plane 
with q = 0.95, Uo = wr, K = 350tJR,iV = 10''. The numbers 
on the plot indicate the number of solutions that exist for 
the steady state photon number in the cavity. Inside the 
swallowtail shaped curve there are five solutions and hence 
multistability, see Fig. |16| 

theory, which we now discuss. As will become clear, the 
structures we see in parameter space are not only generic, 
but are organized in a very particular, and therefore pre- 
dictable, way. 



IX. CATASTROPHE THEORY ANALYSIS 

A. Overview of catastrophe theory 

Catastrophe theory is a branch of bifurcation theory 
that concerns the study of singularities of gradient maps 
[351I37|. Examples of gradient maps abound in physics, 
for example Hamilton's principle of least action in me- 
chanics, and Fermat's principle of least time in optics. 
In both theories the physical paths (rays) are given by 
the stationary points of a generating function $(s;C), 
which in mechanics is the action and in optics is the op- 
tical distance. In both cases the gradient map is obtained 
from a variational principle that takes the form 
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FIG. 16. Plot of multistable steady state photon number riph 
versus rimax (equal to ti^/k^) for — 1630 ojr, and q = 0.95, 
Uo = uj-R., K, = 350 OJR, = 10*. The Ac value is chosen from 
the region which supports five solutions in Fig. |15| 
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FIG. 17. (Color online) Plot of tristable band structure. The 
parameters are given by = QSOujr, Ac — 1640 ojjj, k = 
350 ujR, N = 10*, and Uo = ljr. 
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FIG. 18. (Color online) The cusp catastrophe is generated 
by the quartic potential function <1> given in Table [I] which 
can be viewed as representing a double-well potential. The 
red folded-over surface plotted in this figure obeys the cubic 
state equation d^/ds — s'^ + C2S -I- Ci = 0, which gives the 
stationary solutions s' of $. When C2 < there can be up to 
three stationary points, s^, and s^, for each value of Ci and 
C2. These points are the two minima and single maximum 
of the double-well potential. When C2 > there is only one 
stationary point s corresponding to the minimum of a single 
well. A vertical slice through the figure such that Ci = 
gives a pitchfork bifurcation. The {Ci,C2} plane forms the 
two dimensional control space where the cusp catastrophe it- 
self lives, and this is shown at the bottom of the figure. The 
cusp catastrophe is formed of two fold curves joined at a sin- 
gular cusp point. The cusp catastrophe demarks the region 
of control space that sustains three solutions for s, and so it 
is given by the projection of the folded-over part of the state 
surface onto the {Ci, C2} plane. Crossing the fold lines from 
inside the cusp to outside it, two of the solutions (the maxi- 
mum and one of the minima) annihilate. Mathematically, this 
is described by the potential function being stationary to the 
next higher order gS], i.e. d^^/ds'^ = 3s'^ + C2 = 0. Eliminat- 
ing s by combining this equation with the state equation gives 
the equation for the cusp catastrophe as Ci = ±\/— IGCf /27. 
Right at the cusp point itself, which is given by the control 
space coordinates Ci = C2 = 0, all three stationary points 
coalesce simultaneously to leave a single solution. 



control parameters determine the conditions imposed on 
the paths. For example, if we are interested in the path(s) 
which pass through the point {X, Y, Z} in three dimen- 
sional space, then the coordinates {X, Y, Z} form the 
control parameters. For a fixed set of control parame- 
ters, the potential function defines the height of a land- 
scape with coordinates s. The classical paths (rays) are 
then the stationary points s' (labelled by the index i if 
there are more than one) of this landscape, namely the 
mountain peaks, valley bottoms and saddles [35] . 

The gradient map becomes singular when there is more 
than one stationary point for a given set of control pa- 
rameters. In optics this is the phenomenon of focusing, 
because more than one ray passes through the same phys- 
ical point C = {X, Y, Z}, leading to a caustic. The caus- 
tic, or catastrophe, as it is known in catastrophe theory, 
lives in control space, which in the standard optics case 
is physical 3-dimensional space. As C is varied one can 
explore the structure of the caustic. This is what was 



shown above in Figures 14 and [T5j which are two dimen- 
sional slices through the parameter space {Ac, 77, g, f/o}. 
The crescent and swallowtail shapes are the catastrophes, 
whose full structure only becomes apparent when viewed 
in four-dimensional parameter space. 

Catastrophes are points, lines, surfaces, and hypersur- 
faces in control space across which the number of solu- 
tions to the problem changes. Catastrophe theory clas- 
sifies these catastrophes in terms of their codimension 
K, which is the difference between the dimension of the 
control space and the dimension of the catastrophe. For 
example, if we consider the two dimensional space shown 
in Figures 14 andjlsj we find "fold" curves (K = 1) and 
"cusp" points {K = 2). If we add a third dimension then 
we would find fold surfaces {K ~ 1), cusp edges {K — 2), 
and "swallowtail" points {K — 3). 

In order to make the foregoing discussion more con- 
crete, consider the structure shown in Fig. [18] which il- 
lustrates the cusp catastrophe. The surface shown in 
the figure is the state surface d^/ds — plotted in a 
composite of control and state space. The control space 
is two dimensional and is given by the Ci , C2 plane at 
the bottom of the figure. As listed in Table [ij the cusp 
catastrophe is described by a quartic potential function 
which, by varying the control parameters Ci and C2, 
can be tuned between being a double or a single well 
potential. A prominent physical example of such quar- 
tic potential function is the thermodynamic potential in 
Landau's phenomenological theory for continuous phase 
transitions IMl 



In catastrophe theory, this equation is sometimes referred 
to as the state equation, the generating function <I>(s; C) 
is called the potential Junction and the variables appear- 
ing in the potential function are considered to be of two 
basic types: the state variables s = {si, S2, S3 . . .} and 
the control parameters C = {Ci, C2, C3 . . .}. The state 
variables parameterize all possible paths (not just the 
paths corresponding to the stationary points) and the 



$(s; P, T, h)=^a + Bs* 



As'~ 



hs. 



(35) 



The order parameter s for the phase transition can be 
identified as the state variable. The parameter h de- 
scribes an external field, and A and B are functions of 
pressure P and temperature T. At first sight, it appears 
as though the Landau potential function contains three 
control parameters, and so does not correspond to the 
cusp catastrophe. However, it is easy to see that one of 
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TABLE I. Standard forms of the cuspoid catastrophes 



Name 




K 


Fold 


sV3 + Cs 


1 


Cusp 


sV4 + C2sV2 + Cis 


2 


Swallowtail 


sV5 + CasVS + C2SV2 + Cis 


3 


Butterfly 


sV6 + + Css^/S + C2SV2 + Cis 


4 



the parameters is redundant because the state equation 
can be written in terms of only two control parameters: 
Ci = h/B and C2 = A/B. The Landau thermodynamic 
potential can therefore be seen to correspond to the cusp 
potential function. 

In the present case of ultracold atoms in an optical 
cavity, the potential function is the reduced hamiltonian 
(13), the state variable is the wave function -0, and the 
control parameters are {Ac^rj^q^Uo}. We assume that 
the cavity decay rate k and number of atoms N are 
constants that are unchanged throughout the analysis. 
The stationary Schrodinger equation, obtained from the 
time-dependent Schrodinger equation (10), is therefore 



the state equation which determines the allowed clas- 
sical "paths" (rays) tp. However, because we are in- 
terested here in solutions of the Bloch wave form, our 
"paths" ip are Mathieu functions labelled by the quasi- 
momentum and the depth of the optical lattice. The 
quasi-momentum q is one of the control parameters, but 
the depth of the optical lattice is determined uniquely, via 
the self-consistency equation ([9]), by the number of pho- 
tons in the cavity riph- Therefore, we can choose to work 
with Tiph rather than ^ as already discussed in Section 
|IVB[ The state equation which determines its stationary 
values is the self-consistency equation 

The purpose of formulating our problem in terms of 
catastrophe theory is that we can now take advantage of 
a very powerful theorem |33]. This states that there are 
only a strictly limited number of different forms which 
the potential function $(s; C) can take in the neighbour- 
hood of a catastrophe. The first four are listed in Table|l] 
Note that each of the standard forms is a polynomial in 
s, but is linear in the control parameters. In fact, for con- 
trol spaces of dimension four or less there are only seven 
distinct structurally stable potential functions. Three of 
these require two state variables, whereas we only require 
one state variable, so that leaves the four so-called cus- 
poid catastrophes, which are the ones listed in Table [ij 

This remarkable result allows us to predict, at least 
qualitatively, the structures seen in Figures [14] and [15] 
given only very rudimentary information such as the 
number of control parameters. However, the sting in the 
tail is that it is rare for the potential function that ap- 
pears in any particular problem to already be in one of 
the standard forms shown in Table [l] Rather, it is gener- 
ally necessary to perform various transformations upon 
the variables in order to manipulate the raw potential 
function into one of the standard forms. We already saw 
this for the rather simple case of the Landau theory dis- 



cussed above, and we shall see below that this is also true 
for the problem of multistability in atom-cavity systems. 

From Table ]l] we see that the butterfly potential func- 
tion is the only one which gives up to five stationary so- 
lutions and has a four dimensional control space. We can 
therefore immediately say that our problem corresponds 
at least to a butterfly catastrophe because we have al- 
ready found parameter regimes which give five solutions. 
This does not rule out the possibility of a higher catas- 
trophe (the higher catastrophes contain the lower ones, 
as can be seen in Fig. [18] for the case of the cusp and 
the fold), but until we find regimes with a higher num- 
ber of solutions (and, indeed, we have not) we shall work 
with the hypothesis that we are dealing with a butterfly 
catastrophe. We might therefore expect to flnd a very 
special point in parameter space where all five solutions 
merge into one (the "butterfiy point"). However, this 
requires us to be able to maneuver through parameter 
space in order to find this point. The delicate issue of 
whether the four experimental parameters {Ac, 77,(7, Uq} 
can be transformed into the butterfly's four linearly in- 
dependent control parameters {Ci, C2, C3, C4}, and so 
allow us to fully explore the butterfly catastrophe, will 
be studied in Section IXC below. First, we begin with 



the simpler case of shallow lattices. 



B. Application of catastrophe theory to shallow 
lattices 



The starting point for our analysis will not be the po- 
tential function <I>, which in our case is the energy func- 
tional ( 13 ) that explicitly depends on the atomic wave 
function ■0g(a;, riph) and implicitly on the photon num- 
ber TZph- Instead, we begin one step further along with 
the state equation which, as explained above, is provided 
by the self-consistency equation Eq. ^ for the photon 
number. This was also the approach adopted in [? ] 
in order to tackle bistability in traditional laser systems. 
The photon numbers that satisfy the state equation for 
given values of the control parameters form the set of 
stationary points of the potential function. For our state 
variable we shall actually choose 



V = [/oriph ■ 

In terms of v the state equation can be written 



(36) 



v + v{A,-NUof{v,q)f -V^Uo = 0. (37) 

where (cos^(a;)) = f{v,q) is evaluated in a Bloch state 
ijjqix, v), and the choice of v as the state variable is moti- 
vated by the dependence of the Bloch state on the prod- 
uct C/o^phj which is the depth of the optical lattice, see 
Eq. ([4|. In the above equation we have also rescaled all 
frequencies by k (i.e. we have divided throughout by 
and set k = 1). 

Equation p7[ ) is not straightforward to analyze be- 
cause we do not have a closed-form analytical expression 
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for f{v, q). Thus, we find ourselves in the common situ- 
ation, as mentioned above, that it is not obvious which 
standard potential function, and hence which standard 
state equation, corresponds to our problem. However, 
we have already seen in Sec. |V] that when the depth of 
the optical lattice is small, we can perform a series ex- 
pansion and obtain an approximate analytical expression 
for f{v, q). This is the approach we shall follow first and 
we will take up the general case in the next subsection. 
Specializing to the case of g = 0, we use Eq. ( 26 1 to write 
f(v,q = 0) = 1/2 — w/16, and upon substitution into 
Eq. (|37]) this gives 



+ hiv^ + b2V + &3 = 0, 
32Ac 



(38) 



b2 
bs 



16, 



NUo 

64 (4 + {NUo - 2A,)2) 



-25677^ 



The leading term in the above equation is cubic in the 
state variable v (similar in that respect to the classical 
Kerr non-linearity [21 |4]). It is therefore close to, but 
not yet identical with, the standard form of the state 
equation for a cusp 



+ C2S + Ci = 0. 



(39) 



Complete equivalence to the cusp can be achieved by 
removing the quadratic term from ( 38 1 via the transfor- 
mation s = V + bi/3, leading to the following values for 
C2 and Ci 



- 62 - 

Ci = bs-h,b2 + ^bl 



(40) 
(41) 



Thus, we see that the canonical control parameters in 
the final mapping to the standard form are complicated 
functions of the physical control parameters 77, A^, and 
Uo. 

Let us compare the above prediction to the case illus- 
trated in Fig. [14] The figure shows the number of steady 
state solutions for the photon number in the {77, A^} 
plane, for fixed values of Uo and q. The figure was cal- 
culated for q = 0, and the part of it close to the horizon- 
tal axis corresponds to low photon number, and so the 
shallow lattice theory outlined above applies in that re- 
gion. We indeed find that the first derivative of the state 
equation vanishes at all points along the curves separat- 
ing regions with one and three solutions (this is how the 
curves were computed), which are therefore fold curves, 
while at the point with rj = rjo we find the second deriva- 
tive also vanishes, identifying it as a cusp point where all 
three solutions coalesce into a single solution. We there- 
fore find that catastrophe theory correctly accounts for 
the structure seen in Fig. [Ml 



A key point to note from the above analysis is that the 
underlying catastrophe that we identified had only two 
control parameters even though there were three "exper- 
imental" parameters that could be varied in the original 
statement of the physical problem (we set q = 0). We 
met a similar situation for the Landau theory of continu- 
ous phase transitions discussed above. The question then 
arises, how do we identify the underlying catastrophe in 
cases where the transformation to standard form is hard 
to find? One way to proceed is via the defining character 
of each potential function in Table [T[ which is the highest 
derivative that vanishes at the most singular point. For 
the cusp the most singular point iss = Ci — C2 — 
where the first, second, and third derivatives of the po- 
tential function with respect to s vanish. We shall take 
advantage of this defining character in order to tackle 
the general case of a lattice of arbitrary depth in the 
next subsection. 

It is worth commenting on the role of the extra di- 
mension in control space that was present in the original 
statement of the shallow lattice problem, as given by Eq. 
(38 1. It is easy to imagine that, given a basic catastro- 



phe, we can always embed it in a control space of higher 
dimension without fundamentally changing the catastro- 
phe providing we extend it into the extra dimension in a 
trivial way. For example, given the cusp catastrophe that 
has a minimal control space which is two-dimensional, 
as shown in Fig. |18[ we can always add a third dimen- 
sion such that the cusp becomes a structure rather like 
a tent, with the ridge pole being a cusp edge (a continu- 
ous line of cusps) . The existence of the line of cusps can 
be inferred from the shift of variables s — v + bi/3 that 
was performed: the location of the highest singularity 
is parameterized by the parameter bi which is the extra 
control parameter. 



C. Application of catastrophe theory to lattices of 
arbitrary depth 

Lifting the restriction of small photon number, we de- 
fine the notation Q{v] Ac,r], q,Uo) for the left hand side 
of the state equation ( 37 1 



Giv; A„ r,, q,Uo)=v + v (A, - NUofiv, q)f 
- V^Uo = 0. 



(42) 



The fact that we have four experimental parameters 
holds out the possibility that the above state equation 
is that of a butterfiy catastrophe (see Table |l]) . Further- 
more, because we have already discovered regimes with 
five solutions, we must have at least a butterfiy. How- 
ever, from our experience with the shallow lattice case, 
we know that we may not be able to fully explore the 
catastrophe if some of the control parameters are trivial. 
We shall therefore investigate whether we can find a sin- 
gular point (the butterfly point) where the fifth and all 
lower derivatives of the potential function with respect 
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FIG. 19. A plot showing the values of l/{NUo)^ obtained 
from the simultaneous solution of Eqns 1 43 1-( 45 1 for differ- 
ent values of quasi-momentum. These equ atio ns give the 
first three derivatives of the state equation (42 1, and hence 
correspond to swallowtail points. Only the crosses satisfy 
l/{NUo)'^ > and occur only for q > 
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FIG. 20. (Color online) A plot showing the values of the state 
variable v — Uoriph for which the first three derivatives of the 
state function Q vanish simultaneously (swallowtail points). 
For approximately 0.6 < q < 0.8 there are two such points for 
any given q. 



to V simultaneously vanish (i.e. the fourth and lower or- 
der derivatives of G{v; Ac,ri,q,Uo) simultaneously van- 
ish) and Eq. (42) is also satisfied. 



Taking the derivatives of Eq. ( 42 ) , and simplifying 



slightly the resulting equations, we find 



\NUo \NUo 



-/ = 



2f + vf" NUo 

3//" + 3(r)2 + V (arr + fn) _ Ac 



+ NUo 

yjiv _|_ 4Jiii 



(43) 
(44) 

(45) 

A, 

NUo 
(46) 



where the first, second, third, and fourth derivatives of 
the function f{v,q) with respect to v are denoted by 
respectively. In fact, we have seen Eq. 



combined into a single equation: 

2/r + t'(/r + (/')') 

3//" + 3(r)2 + z;(3/v" + /r") 







(47) 



We solve this equation for v at different values of q by 
numerically finding the zeros of the left hand side. Once 
we find the zeros for a particular g, we can use Eq. (44 1 to 



calculate Ac/NUq at these values. Next, these values of v 
and Ac/NUo are used in Eq. ([43]) to calculate l/{NUo)'^. 



In this last step we find that we obtain the unphysical 
result l/{NUoY < unless q > qsw, where ^sw = 0.545 
is a certain critical value of the quasi-momentum. This 
is illustrated in Fig. [19] where we plot the values of 
l/(-/V[/o)^ computed using the above method for differ- 
ent values of the quasi-momentum q in the neighborhood 



of 



Only the crosses satisfy l/{NUo) > 0. We 



drop the solutions corresponding to the dots, for which 
1/{NUq)'^ < 0. For q > qsw there are always values of 
q at which l/(A^C/o)^ > 0. Note that g^w = 0.545 is a 
universal result since it does not depend on any other 
parameter values. 



(43) before, as it is the same as Eq. (25) that we used as 



The final step is to check if Eq. ( 46 ) for the fifth deriva- 
tive is satisfied at the values of Ac, NUq, and v that we 
computed using the lower derivatives. We did not find 
any value of g > gsw where this was the case. The im- 
cilitate this, observe that Eq. (|44|) and Eq. (|45|) can be portant part of our numerical computation involves the 



a condition for bistability. Our strategy will be to find 
solutions to Eqns ( 43 - 46 ) numerically. In order to fa- 
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calculation of the derivative of the function f{v,q) for 
which we do not have an analytical expression. We used 
the MATLAB® gS] routine DERIVSUITE [SD] to cal- 
culate the derivatives. This routine also provides errors 
on the derivatives and we can compound errors and find 
values for expressions like the left hand side of Eq. ( 47 ) 
with error. We use this error as the tolerance in our zero 
finding. Some more details of this procedure are provided 
in the Appendix. 

The fact that we did not find a point where the four 
higher derivatives of the function Q simultaneously van- 
ish in the range < g < 1 (we do not consider negative q 
because fiv,q) is symmetric under g — >■ — g) means that 
although the underlying catastrophe that organizes the 
solutions is at least a butterfiy (because we find five solu- 
tions), the four experimental parameters at our disposal 
{Ac, ?7, q, Uq} do not translate into four linearly indepen- 
dent coordinates in control space {Ci, C2, C3, C4} [5T]. 
We are therefore not able to navigate freely through the 
four-dimensional control space and locate the butterfly 
point at Ci = C2 = C3 = C4 = 0. This is the extension 
into four dimensions of the situation we already found in 
Section liXB I for shallow lattices. 

The identification of places where three derivatives of 
G simultaneously vanish means that the highest singu- 
larities we have in the parameter space {A^, rj, q, Uq} are 
swallowtail points (swallowtails are contained within a 
greater butterfly catastrophe). In the Appendix we out- 
line a proof that in the neighbourhood of a point where 
three derivatives of the state equation vanish the poten- 
tial function must be equivalent to that of a swallowtail 
catastrophe. Note that these swallowtails occur entirely 
in control space and are thus true swallow tail catas- 
trophes in the sense of Table IT] The swallowtails shown 
previously in Figures |3j [9j [TOj [H] [12] and [13] are not swal- 
lowtail catastrophes because these figures show a com- 
bination of state and control spaces. By contrast. Fig- 
ures [MJ [15] and [21] are solely in control space. In par- 



ticular. Fig. 21 shows a two-dimensional slice through 
the three-dimensional swallowtail catastrophe in control 
space. Unlike Fig. [15] this slice includes the swallowtail 
point, which is the highest singularity on the swallowtail 
catastrophe, and is the point where four solutions of the 



state equation ( 42 ) simultaneously coalesce so that num- 
ber of solutions changes by 4 (i.e. the point where three 
derivatives of G simultaneously vanish). We emphasize 
that this can only occur when q > qsw We also note from 
Fig.[20]that between 0.6 < q < 0.8 we find more than one 
swallowtail point for a given q, which is again an indica- 
tion of the presence of a higher underlying catastrophe. 
As described in reference [52], swallowtails contain two 
cusps, and butterflies contain two swallowtails. 



X. STABILITY ANALYSIS 
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FIG. 21. (Color online) Bifurcation structure of the solutions 
to the self-consistency equation Eq. ([9| in the {?7, Ac} plane 
with q = 0.96, Uo = 1.13cjr, k = 350tJR,iV = lO"*, and 
hence NUo/2 = 16.1k. The numbers on the plot indicate the 
number of solutions for riph in that region of the parameter 
space. The inset shows a swallowtail singularity point where 
five solutions coalesce into a single solution. The coordinates 
of the swallowtail point shown in this figure are v — 0.04, 
77 = 1.7k, Ac = 6.4fi:. 



approach more in line with 
tional, Eq. 



H], where the energy func- 



(13), and the nonlinear equation of motion, 
Eq. (10 1, for the atomic wave function serve as the start- 
ing points for an examination of energetic, and dynamic 
stability, respectively. Hence, we examine the stability of 
the Bloch states at different values of quasi-momentum 
at fixed values of f] and Ac. Before going into the details 
of the calculation we note that it is well known from the 
study of bistability in classical nonlinear cavity systems, 
that the back-bent branch of the lineshape profile shown 
in Fig.[2]is unstable, see, for example, |3]. From the three 
dimensional plots Fig. 12 and Fig. [13] we see that this 
back-bent branch corresponds to the upper branch of the 
loop in energy-quasi-momentum space. Thus, we expect 
to find this branch unstable. 



The stability of cold atoms in an optical cavity has 
been treated in Refs. [TUISni. In this case we follow an 



We first consider energetic stability in the spirit of [TH] ■ 
The grand canonical potential per unit length is G[tp] = 
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cm 1 




NUo l-T^ 



/J'dxiV'l cos2(a;) 



/i / dx\'ip{x)\' 



(48) 



We perturb the wavefunction as "ipi^) — 4'q{x) + Sip{x), 
where ipo extremizes G, i.e. one of the solutions that 
we obtained in Sec. |V] Since ipQ is an extremum, the 
first order variation in G vanishes and the second order 
contribution can be written as 



N 



+ 2p{Si;\ cos2(z)|i^o)(i^o| cos^ {z)\6i;) , 



where 



and 



Ho 



+ [/o^ph cos^(z) 



A,-NUo(iPo\cos'^{z)\iPo) 



Ae-AfC/o(V-o|cos^(z)|V.o) 



(49) 



(50) 



(51) 



Equation ( 49 ) can be cast into a simple matrix form 

dx'^\x)A^{x). 



6G2 ^ 1 
iV ~ 2 



(52) 



Here 



and 



A 



Ha+2pcos^(x)i,o(x)I*[.. 
2pcos^(x)^l{x)r[...] 



2/ocos {x)'tl>o{x)I[. . .\ 
ffo+2pcos^(a;)V>o(a;)7[. 



where /[.••] is an integral operator defined by 

I[5^*{x)] = ( dxcos^{x)%{x)5^j*{x). (53) 



The eigenvalues of the matrix A decide the energetic sta- 
bility. If A is positive definite, the solution i/jo is energet- 
ically stable. 

In order to examine dynamic stability, we linearize the 



equation of motion Eq. (10) by writing ?/;(a;, t) = [il^oix)- 



5'^|}{x,t)]e-'■^'\ This leads to 



Ax' 



+ ?7o I ao P cos^(a::) - p 



5iIj{x) (54) 



■ 2/ocos2(x)V'o(a;) (r[(5i/>(x)] -f I[5i1j*{x)]) . 




-1 -0.5 0.5 1 

quasi-momentum q (units of kc) 

FIG. 22. (Color online) Energetic and dynamical stability of 
the band structure loops. The upper branch of the loop (black 
dashed line) is energetically and dynamically unstable. The 
other branches (red lines) are energetically and dynamically 
stable. Parameters are, Ac — 3140 a;R, k = 350 o^r, Uo = (jJr, 
N = 10", 7? = 2.8 ryo, and -qo = 325 ojr. 



One can write a similar equation for Stp* and combine 
the two into a matrix equation similar to Eq. ( 52 1 



at 



(55) 



where cr^ is the Pauli z-matrix. The solution ipg is dy- 
namically stable if all the eigenvalues of (TzA are real. 
Thus, the occurrence of complex eigenvalues of azA sig- 
nals dynamical instability. Before we quote the results, 
a comment is in order about the form of the pertur- 
bations Sijj. The integral operator in Eq. (53 1 cou- 



ples the perturbation and i/jq. If i/jq = e^'^^Uq{x), with 
Uq(x) = Uq{x + tt), i.e. a Bloch function with quasi- 
momentum q, the form of Sijj that leads to non-zero cou- 
pling is 



5'ip{x) 



3 



^i2jx 



(56) 



That is, the perturbation should be a Bloch wave with the 
same quasi-momentum as ipQ. In Eq. (49 1 we consider the 



change in the grand canonical potential per unit length, 
but the above choice is made to satisfy the requirement 
that the integral operator / over the system size gives a 
non-zero answer. A physical way to motivate the above 
choice goes as follows: in the absence of interatomic inter- 
actions the allowed excitations of the quasi-momentum 
state q have to be in multiples of the crystal momen- 
tum 2kc (simply 2 in dimensionless units) since the only 



23 



source of perturbation is the interaction with the cavity 
field. The above form respects this requirement. Also, 
the number of terms in the Fourier expansion Eq. ( 56 ) 



has to be less than the number of terms in the original 
expansion for tpo in Eq. ( 16 ) to avoid spurious instabili- 
ties [18] . Using Eq. ( 56 ) , we find that the upper branch 



of the looped dispersions is both dynamically and ener- 
getically unstable as expected. The other two branches 
are stable. This is shown in Fig. 22 for one particular 



case. We shall not perform the stability analysis for the 
case when there are five solutions, but anticipate by an 
extension of the case for the bistability scenario, that two 
of the solutions will be dynamically unstable and three 
will be stable. 



XI. SUMMARY AND CONCLUSIONS 

In this paper we have analyzed bistability in atom- 
cavity systems in situations where the atoms are in spa- 
tially extended states (Bloch waves) with non-zero quasi- 
momentum q. We find that bistability in the number of 
photons in the cavity goes hand-in-hand with the emer- 
gence of loops in the band structure. Both are manifes- 
tations of a bifurcation in the stationary solutions of the 
coupled atom-light equations of motion. 

We have studied how the loops appear and disappear 
as the laser detuning and the laser pumping rate are 
changed. In particular, Eq. (31) provides an analytical 



estimate of the critical pump strength 77cr(<?) at which 
bistability sets in. It depends on the quasi-momentum of 
the atomic state, and predicts that loops first appear at 
the edges of the first Brillouin zone {q = ±1) and then 
move inwards. This is indeed what we find upon solving 
the coupled atom-light equations numerically: swallow- 
tail loops appear at the edges of the first Brillouin zone 
as the pump strength 77 is increased above ?7cr(9 = !)• As 
77 is increased further the swallowtails extend inwards, 
merge, and detach from the rest of the band to form a 
separate loop centred at g = which ultimately closes 
up and vanishes. A rather similar behaviour is observed 
as the pump-cavity detuning Ac is swept from below the 
cavity resonance to above it. 

The loops we find are qualitatively different from those 
that occur for BECs in static optical lattices in the pres- 
ence of repulsive interatomic interactions [T7HI1]. There, 
the loops are centered at the edge of Brillouin zone and 
cause the dispersion to have a finite slope at that point. 
By contrast, the band structure we find always has zero 
slope at the edge of the Brillouin zone. Nevertheless, 
there are also many similarities, including the stability 
of the various branches of the loops. We find that the 
upper branch of the loops are energetically and dynam- 
ically unstable, as expected from optical bistability con- 
siderations. 

The extra degree of freedom afforded by the quasi- 
momentum (over considering only q = 0) results in the 
possibility of tristability, namely regions of parameter 



space where there are five solutions, three stable and two 
unstable. The complexity of the solutions in parame- 
ter space led us to perform an analysis of the problem 
in terms of catastrophe theory which is a useful mathe- 
matical tool for understanding the organization of bifur- 
cations of solutions. The key to our treatment was the 
recognition that, because exact solutions for the atomic 
wave functions are Mathieu functions which are speci- 
fied only by the lattice depth (once one chooses a quasi- 
momentum), the photon number Uph which determines 
the cavity depth provides a completely equivalent de- 
scription to the wave function. 

In the case of shallow lattices we were able to proceed 
analytically and found that the structure of the solutions 
in parameter space when q — corresponds to a cusp 
catastrophe, at the most singular point of which three 
solutions (two stable and one unstable) form a pitch- 
fork bifurcation, and this describes the onset of bista- 
bility as the laser pumping is increased. Interestingly, 
the three experimental parameters {Ac,?7, C/o} reduced 
to just two effective control parameters. In the general 
case of arbitrary lattice depth and < g < 1, the high- 
est singularities we found were swallowtail catastrophes 
where four solutions simultaneously merge. The swallow- 
tails only exist when q > qsw — 0.545. However, there 
is good evidence that there is an underlying butterfly 
catastrophe, but, once again, the experimental parame- 
ters {Ac,77, g, C/q} reduced to three effective control pa- 
rameters meaning that generically one is unable to locate 
the butterfly point (where five solutions simultaneously 
merge) . 

The band structure loops found here have important 
implications for Bloch oscillations of atoms in cavities 
[5U1 [5T] . Bloch oscillations are essentially an adiabatic ef- 
fect where, as a force is applied to the atoms they remain 
in the same band but their quasi-momentum evolves lin- 



early in time, as shown in Eq. (22). Swallowtail loops in 



the band structure will have a deleterious effect on Bloch 
oscillations because, as the quasi-momentum evolves, the 
atoms will reach the edge of a loop where the branch 
they are following vanishes. This will lead to a sudden, 
non-adiabatic, disruption in the state of the atoms as 
they are forced to jump to another branch or even an- 
other band. For BECs in ordinary static optical lattices 
these non-adiabatic jumps are thought to be the cause of 
the destruction of superfluidity during Bloch oscillations 
(52 . We have not included interactions in our treat- 
ment (interactions are necessary for superfluidity), but 
related effects will likely occur, especially considering the 
added heating due to the fact that the lattice depth will 
also abruptly change at the same point. However, when 
the loop detaches from the main band it will no longer af- 
fect Bloch oscillations. Furthermore, loops only occur for 
certain limited regions of parameter space, i.e. inside the 
cusp and swallowtail catastrophes shown in Figs Mj 15 
and |21[ For experiments involving Bloch oscillations we 
therefore recommend that parameter regimes are chosen 
which lie outside to these regions. 
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Finally, we add that although we have only considered 
Bloch waves in this paper, localized states (for exam- 
ple Wannier functions) can be formed from superposi- 
tions of Bloch waves with different values of the quasi- 
momentum. In this sense, localized states therefore con- 
tain all values of the quasi-momentum and so might be 
expected to display tristability too. However, it should 
be borne in mind that the nonlinearity of the system 
means that superpositions of Bloch states of different q 
but the same lattice depth will not in general obey the 
effective Schrodinger equation (10). 



ACKNOWLEDGMENTS 

We gratefully acknowledge E. A. Hinds, D. Peli- 
novsky and M. Trupke for discussions. For funding, 
DHJO and BPV thank the Natural Sciences and Engi- 
neering Research Council of Canada, and the Ontario 
Ministry of Research and Innovation, and JL thanks 
VR / Vet enskapsradet . 



Taylor expansion is of order fc -I- 1), can be locally 
expressed as p{y{s)) with y: 7?." — )■ TZ^ being a 
smooth reversible change of co-ordinates. 

• is the vector space of polynomials in si, . . . , s„ 
of degree < k. 

• J^l is the subspace of with zero constant term 

• Afc(p) is the subspace of J,^' spanned by all 

Qf [§rj ' where 1 < i < n, Q e E'^^, and the 
bar symbol represents the restriction to the fc-th 
order of the expansion. 

• The codimension of a function p is the codimen- 
sion of Afe(p) in for any k for which p is k- 
determinate. 

• An r — unfolding of p at is a function: 



P : 7^"+'■ ^ 7^, 

(Sl , . . . , Sji ; ^1 ; • 



APPENDIX 

We shall now sketch a proof showing that the function 
Q defined in Eq. ( 42 ) produces swallowtail catastrophes 



between 0.545 < g < 1 (in fact it produces two lines of 
swallowtail points, as shown in Fig. 20 1 where its deriva- 



tives up to third order vanish simultaneously |361 154j . 

Consider, for example, the point in control and state 
space given by Cq = {A^ = 0.90, rj = 14.5, q = 0.69, Uq = 
0.15} and wq = 7.75, riph = 50.3 (frequencies are mea- 
sured in units of k and the number of atoms is set at 
N = 10^). The numerical package |SD] we used for cal- 
culating the derivatives gives error bounds allowing us to 
estimate the accuracy of our calculations. For the point 
{Co , Wo} we found that the right hand side of Eq. ( 47 ) was 
equal to 8.09 x 10"^^ with an error of 4.00 x lO^^This 
means that the third derivative of the state function van- 
ishes within error, indicating a swallowtail point. How- 
ever, the smallest value we found anywhere in parameter 
space for the difference between the left hand side and 
the right hand side of Eq. (461 was —8.64 x 10~^ with 
error 6.00 x lO^^'^. This means that the fourth deriva- 
tive of the state function does not vanish within error, 
suggesting there is no butterfly point. The value of the 
quasi-momentum at the point where we found this min- 
imal fourth derivative was q — = 0.545. 

Before outlining the proof, we first give some defini- 
tions |36j . If n is the number of state variables, then 
consider a function p: 7?." — >■ TZ 

• j^p is the Taylor expansion of p to order k 

• J^p is i^'p minus its constant term 

• p is fc-deteminate at if any smooth function p + q^ 
where q is of order A: -|- 1 (leading order term of the 



such that Po,...o(s) —p{s)- 

More informally, the term "unfolding" refers to how the 
catastrophe unfolds as one moves away from the origin 
in control space. At the origin in control space the catas- 
trophe reduces to its most singular part, known as its 
germ. For example, from Table IT] the germ of the swal- 
lowtail catastrophe is given by s . The terms in the po- 
tential function which depend on the control parameters 
are called the unfolding terms, and the number of them 
is equal to the codimension. If P is an r — unfolding of 
p, set 



H{p), 



.,w, 



dt, 



(j'^(Po,...,t„)) . (57) 



W''{P) is the subspace of J,j spanned by 
{w'liP),...,wHP)}. 

Referring to the standard forms given in Table [T) the 
potential function and state equation for the swallowtail 
are given by 



Hs;C) = -s' 
5 



C: 



Q{s;C) 



+ Css^ + C2S + Ci 



Cis 
= . 



(58) 
(59) 



Notice that the state equation for a swallowtail catastro- 
phe is similar to the potential function for a cusp catas- 
trophe up to a constant Ci. Since the state function Q 
is the central object in the treatment given in Section 
|IX| rather than the potential function, instead of proving 
that the underlying potential function is equivalent to 
that of a swallowtail catastrophe, we will prove that the 
state function Q around the singular point vq and Cq is 
equivalent to the potential function of a cusp catastrophe 
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(note that this is different from the small photon number 
case we studied in Section IIXBI where we showed that 
the underlying potential was equivalent to the potential 
for a cusp catastrophe). To that end, first notice that 
the role of the constant term Ci in Eq. ( 59 ) is played by 
—ti^Uq in Eq. (42 1. Subtracting this function we have 



a modified form of Q (where we have also dropped the 
dependence on q since we are focusing on a particular 
quasi- momentum) : 



Fiv; {A„ Uo}) = F{v; C) ^ v + i;(A, 



NUof{v,q)f 
(60) 



which satisfies F {vo;Cq) — F {vq;Co) — F (vq^Co) — 
0. In order to have a function defined in the neighbor- 
hood of vq and Cq, let us set the origin of v at Vq and the 
origin of C at Co and define 

F,{v; {A„ Uo}) = Fiv + vo, C + C'o) - F{vo, Cq). (61) 

Thus, we have i^i(0,0) = and for the function g{v) = 
Fi{v, {0, 0, 0}) the most singular point is at w = where 
g , g , g vanish. The function g is the germ which we 
described above, and is the key feature which identifies 
the catastrophe. When g is Taylor expanded around 
one has 



9iv) = -^v 



(62) 



where ^'•"•'(O) is the first non-zero Taylor coefficient. This 
means that g is 4-determined around and we say that 
g ^ around 0. According to Table |lj the canonical 
unfolding of the 4-determined germ around is the cusp 
catastrophe $(s; C) — s^/4 -I- C2S^/2 -|- Cis, where v and 
s are related via a differomorphism (smooth transforma- 
tion of coordinates). 

Next we calculate the codimension of g. The Jacobian 
ideal for g is in this case A4{g) — {v'^,v^ + j^r^; l^^o^^^l- 
Hence, the codimension of g is dim{jf) — dim(A4((7)) = 
4—2 = 2. The function Fi is thus a 2 parameter unfolding 
of the germ g. In order to prove that the function Fi can 
be described by a cusp catastrophe, we need to prove 
that Fi is isomorphic as an unfolding to the canonical 



form $(s; C) = + C2s'^/2 + Cis. In order to do this 
we need to invoke the idea of transversality. 

Transversality generalizes what we know of two inter- 
secting lines in a two dimensional plane to multidimen- 
sional manifolds. Two subspaces of a manifold are trans- 
verse if they meet in a subspace that is as small in di- 
mension as possible. If Xi (dim r) and X2 (dim t) are 
subspaces of X (dim n), Xi and X2 are transverse if 
their intersection is empty or if it is of the dimension 
max(0,r + t — n). Our first aim is to prove that the 
2 — unfolding of Fi is a versal unfolding. To do this we 
use a defining theorem for versality from |36j which states 
that: an r — unfolding P of p, where p is fc-determinate 
is versal if and only if W''{P) and Afe(p) (defined above) 
are transverse subspaces of J^. We have already found 
A4{g), the polynomial space W^^(Fi) is spanned by the 
vectors: 



Wi(F,) = 

W2{Fi) = 



^^{JHF,{0,Uo,Q,))), 
^(j^Fi(0,0,A,))), 



The expressions depend on the derivatives of the coupling 
function /(u, q) and the value of the parameters at the 
singular point {vo,Co}. They are too cumbersome to 
state here but their general forms are given by 



which we determined numerically and all of the z^^'s 
are non-zero. The polynomials Wi are linearly indepen- 
dent which gives the dimensionality dim(Vt^''(i^i)) = 2. 
Furthermore, we have verified that the rank of the ma- 
trix formed by the polynomial coefficicints of A^^g) 
and W^^(i^i) is 4 and this combined with the fact that 
dim(A4(5)) + dimiW^iFi)) ^2 + 2 = dim( Jj*) proves 
that A4((7) and W'^iFi) are transverse. Thus, by the the- 
orem stated above F\ is a versal unfolding of the germ g 
and since it is a 2 — unfolding (codimension of 5 = 2) it 
is also universal [36] . This proves the equivalence of the 
unfolding of F^ to the cusp catastrophe. 
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